Tanager Evolution Project Notes
2024-08-13
initialization
suppressMessages({
ddPath <- "/Users/michaelalfaro/Dropbox/color_evolution_scratch/tanager_comp_analysis/allison_data/RefDataSpeciesSummary_MatchingMF"
dat <- pblapply(list.files(ddPath, full=T), read.csv)
wl <- dat[[1]][, 1]
dat <- sapply(dat, "[", -1)
specs <- as.rspec(cbind(wl, do.call(cbind, dat)), lim=c(300, 700))
specs <- procspec(specs, fixneg="zero")
male_specs <- specs %>% select(!ends_with("f"))
crown_data <- male_specs %>%
select(wl, starts_with("Crown_")) %>%
rename_with(~ gsub("m$", "", .)) %>%
rename_with(~ gsub("^Crown_", "", .))
rgb_colors <- spec2rgb(crown_data)
names(rgb_colors) <- names(crown_data)[-1]
head(crown_data)
})## wl AcaBai AniIgn AniLac AniMel AniNot AniSom BanArc BanAur
## 1 300 2.356690 1.1447721 11.39746 8.896685 4.535676 2.014334 5.533874 3.126775
## 2 301 1.903805 0.8218703 12.11596 8.899757 5.386252 1.444982 4.656108 3.745577
## 3 302 1.777677 0.3635018 12.44752 9.508270 6.579554 1.646972 4.473581 4.624964
## 4 303 1.956667 0.5506667 12.46547 9.872667 6.587833 1.983947 4.986000 4.196333
## 5 304 2.086104 0.5410315 12.57834 9.691481 6.914120 2.683912 5.666046 5.145111
## 6 305 2.553706 0.9282766 12.29238 10.253288 7.446329 2.923659 5.490410 6.381721
## BanEdw BanMel BanRot ButAur ButExi ButMon ButWet CalCoc
## 1 1.4074342 2.792863 2.720081 12.36975 9.853977 0 1.940090 0
## 2 0.8987856 2.613824 2.721775 12.86508 9.873162 0 2.141559 0
## 3 1.0957351 2.701293 2.913748 13.68989 10.111072 0 2.995856 0
## 4 1.6188667 2.551333 2.243000 13.83140 10.321933 0 2.810667 0
## 5 2.0386333 2.740259 1.734444 13.45717 10.115626 0 3.157426 0
## 6 2.5069297 3.211658 1.078811 14.04848 10.483344 0 3.042622 0
## CamHel CamPal CamPar CamPau CamPsi CatAna CatDia CatHom
## 1 3.469025 6.643459 2.719078 5.499726 2.579595 4.518847 2.219834 5.587705
## 2 3.036942 7.602843 3.744938 5.698823 3.010506 4.336070 2.382877 6.276141
## 3 2.725917 7.989854 4.219054 5.927595 3.078775 4.244808 2.497530 6.854827
## 4 2.889133 7.241067 3.745280 5.383200 3.367667 3.862467 1.883267 6.427067
## 5 3.143063 6.617119 3.128248 5.333707 3.731115 4.145883 2.269567 5.481470
## 6 3.869150 6.013214 3.264869 5.477861 3.516852 4.021309 1.918708 6.225052
## CatIno CerOli ChaEuc ChlCal ChlNit ChlPho ChlRie ChlSpi
## 1 3.008490 1.597151 5.387928 2.630297 0 1.787583 2.784198 1.612941
## 2 2.616133 1.730229 6.023412 3.474346 0 2.044171 3.276139 1.277324
## 3 3.294292 2.450377 6.458176 3.473768 0 1.696432 3.689607 1.380377
## 4 4.443133 2.863067 6.054750 3.991000 0 2.461333 4.117933 1.372000
## 5 3.437819 2.587056 5.016699 3.586952 0 3.661704 4.370585 1.306522
## 6 5.182589 2.042841 5.806554 3.673250 0 3.939279 4.340295 1.892305
## ChrChr ChrSal CisLev CneRub CoeFla ComBae ComGar ComLor
## 1 1.949883 4.471500 6.664564 4.975148 3.760420 4.967959 1.844345 4.560399
## 2 2.329393 4.400425 6.275407 4.706490 4.126041 5.326747 2.188339 4.799556
## 3 2.523396 4.784778 6.921449 5.457355 4.167218 4.105865 1.724850 5.719318
## 4 2.792556 4.523722 6.041067 6.047067 4.111000 3.140896 2.095556 5.047583
## 5 2.816093 3.580108 6.493004 6.208793 4.068874 3.156895 2.492580 4.814028
## 6 2.818090 3.429653 6.855703 6.840187 4.819775 2.791245 2.552961 4.056009
## ConAlb ConBic ConCin ConFer ConLeu ConMar ConRuf ConSit
## 1 10.46993 2.015616 3.804209 1.9168577 2.578745 4.845310 3.141086 2.525085
## 2 10.12310 1.258997 3.714779 1.8302739 2.849072 4.560984 2.548177 2.491962
## 3 10.40977 1.210195 3.284392 0.9475532 3.220024 4.541942 2.517834 2.725701
## 4 11.18733 1.583111 3.498333 0.9898000 2.851333 4.683467 2.070467 2.394800
## 5 11.26497 1.388296 4.440231 1.8609481 2.321988 4.088907 2.744841 2.564296
## 6 12.11861 1.712550 4.598439 1.9402180 1.952958 3.663535 2.913573 2.412890
## ConSpa ConSpm ConTam CorCuc CorMel CorPil CreDen CreVer
## 1 3.358489 0 3.412506 6.490117 4.495232 10.35447 2.840637 4.538499
## 2 3.709069 0 3.201041 7.144863 4.478360 10.73900 2.686527 4.503564
## 3 3.925063 0 3.037099 7.932816 4.588965 10.93741 1.949874 4.403225
## 4 4.112333 0 3.636267 7.416467 4.771722 12.02407 1.951083 5.057800
## 5 4.758173 0 4.172085 7.681907 5.639179 10.85867 1.922611 4.922030
## 6 4.390084 0 4.023079 7.862213 5.623269 10.67903 1.630532 5.197043
## CyaCae CyaCys CyaLuc CyaNit CyiCyi CypHir DacAlb DacBer
## 1 47.82842 41.49798 23.52521 36.24909 9.682599 3.422173 30.74110 10.52712
## 2 48.71184 42.43106 24.20018 37.04390 10.122437 3.426750 30.94648 11.42346
## 3 49.29641 43.49012 24.71517 37.55738 10.289468 3.745732 31.15261 10.78989
## 4 49.57293 44.47547 24.95913 37.89383 10.688333 3.887600 31.72800 11.02356
## 5 50.42770 45.24845 25.83192 39.48137 10.758269 3.550522 32.22233 11.77883
## 6 50.60203 45.97401 26.38695 40.60357 10.539698 3.741910 32.99691 11.63292
## DacCay DacFla DacLin DacNig DacVen DacVig DelCas DigAlb
## 1 16.62614 7.034602 22.20075 20.04237 24.05480 20.14063 12.71208 5.808184
## 2 17.32347 7.183602 22.53657 20.95381 24.02077 20.86264 13.64604 5.711144
## 3 17.80806 8.018494 23.07261 20.84107 23.55332 21.15301 13.70460 5.501305
## 4 18.50593 8.506333 24.06907 20.39833 23.70200 21.38573 13.86287 5.799800
## 5 19.93100 9.992485 24.85467 20.41631 24.27682 21.83977 14.96010 6.483400
## 6 19.51759 10.959667 26.11854 21.03806 24.74464 22.36942 15.42640 6.707083
## DigBar DigBru DigCae DigCar DigCya DigDui DigGla DigGlo
## 1 6.442543 5.176213 4.761571 3.826953 10.33006 7.791004 12.56828 0.00000000
## 2 5.508182 5.402843 5.634023 3.892699 10.60421 8.172086 12.80130 0.00000000
## 3 5.339232 5.436029 5.975823 2.957868 10.81254 7.808751 12.78159 0.03230631
## 4 4.667250 5.552800 6.400200 3.486667 11.47007 6.575733 13.07420 0.00000000
## 5 5.529088 5.609878 6.701107 4.228474 11.73469 6.108693 13.32874 0.41992963
## 6 5.729264 5.660191 6.361795 3.445634 12.32895 5.910584 13.70849 0.95878018
## DigGls DigHum DigInd DigLaf DigMaj DigMys DigPlu DigSit
## 1 5.911885 5.315526 3.912009 4.740884 5.746431 6.459852 9.968993 7.514505
## 2 5.905845 5.477690 3.296104 5.007212 5.998514 5.987995 10.236485 8.075579
## 3 6.023903 5.611667 3.129428 4.902143 6.072998 6.033641 11.279784 7.875029
## 4 6.213250 5.820533 2.860333 5.006771 6.509533 6.716333 11.134667 8.076500
## 5 6.421384 6.156630 2.969861 4.906539 5.850556 6.852841 11.711237 7.907699
## 6 6.629507 6.380209 3.943491 4.563189 5.928620 7.607613 12.199932 7.886867
## DigVen DiuDiu DiuSpe DolFri DonAlb DubTae EmbDui EmbHer
## 1 4.413730 5.754768 4.513088 5.628032 4.645562 9.388712 1.686559 4.040586
## 2 4.292516 6.051697 4.685247 6.604378 5.038362 9.714932 2.270955 3.951989
## 3 4.044523 5.380742 4.342351 7.106216 4.673121 10.139378 1.531207 3.911493
## 4 3.733833 4.769333 4.060133 6.991333 4.297933 10.404267 3.147333 4.198333
## 5 4.270157 3.463952 4.665000 6.549444 3.256133 10.792793 3.227259 3.965176
## 6 4.472887 2.724789 4.475528 5.723248 2.586216 10.657593 3.542144 4.317468
## EmbLon EmbPla EmbYpi EucPen EunCam GeoCon GeoDif GeoFor
## 1 2.338532 4.269550 4.909153 0.4650937 8.777310 4.302198 4.160865 5.656685
## 2 3.439315 3.636867 4.307609 0.8571243 9.311593 3.814930 4.490456 5.751912
## 3 3.362640 4.218802 3.409551 1.2577694 10.058741 4.273697 4.224193 6.028177
## 4 3.698667 3.939067 4.183800 1.4741333 10.300400 3.505800 4.474000 6.035600
## 5 3.461648 3.932163 5.405856 1.1612111 11.166744 3.235274 4.708181 6.405648
## 6 3.409676 4.654184 6.260063 0.9784072 10.824191 3.115571 4.505560 6.617688
## GeoFul GeoMag GeoSca GubCri HapRus HapUni HemAtr HemCal
## 1 6.448760 6.531546 4.561547 7.380761 4.582564 7.463453 2.444315 1.789450
## 2 5.698004 7.064521 4.258554 6.909324 5.085793 7.039780 2.618054 1.943099
## 3 5.313429 7.979640 4.555024 6.311311 5.553897 8.032144 2.213946 2.463144
## 4 5.100000 8.089267 3.589333 6.047167 5.761733 8.117767 2.134833 2.643067
## 5 4.763537 8.176015 3.280883 5.714065 5.783456 8.154983 2.177606 3.226993
## 6 5.055092 7.564323 2.848854 6.350489 5.951241 7.519610 1.680514 3.041022
## HemFla HemFro HemGoe HemGui HemMel HemPar HemPiu HemRey HemRuc
## 1 3.287751 0 5.605915 3.379796 3.736193 3.1149820 0 1.266595 4.524636
## 2 3.218171 0 5.699159 3.965413 3.895222 3.4821216 0 1.997243 4.850843
## 3 2.946195 0 4.831592 4.463416 3.758582 1.4672354 0 2.621126 5.460076
## 4 2.740067 0 4.679933 4.512533 3.976467 0.9297917 0 3.241667 5.154133
## 5 2.694589 0 4.609824 4.911848 4.363937 0.8993449 0 3.903537 3.785207
## 6 2.550214 0 4.942499 3.846759 4.760578 1.1250101 0 4.106505 3.929441
## HemRus HemSup HemTri HemVer HemXan HetRub HetXan IdiBra
## 1 2.502548 2.098243 3.904831 3.234018 3.045542 6.152712 4.228811 3.312871
## 2 2.309180 1.753862 4.316602 3.620662 2.942859 6.544171 4.682160 3.618069
## 3 2.877843 2.183583 5.085332 4.342088 3.038568 6.654640 4.674836 4.216775
## 4 3.275933 2.895000 5.113800 4.343833 3.104200 6.115667 4.887733 4.643333
## 5 3.262852 2.984352 5.068830 4.830097 3.556463 6.642037 4.921104 4.764395
## 6 3.059625 3.220811 5.025465 5.040315 3.916326 6.796928 4.661458 4.698480
## IncLae IncOrt IncPer IncPul IncWat IriAna IriJel IriPor
## 1 4.414038 5.292791 2.551869 3.270497 4.785303 5.771081 3.992751 4.377826
## 2 3.989401 6.119745 2.390995 3.248391 5.037771 5.888508 4.252169 3.813357
## 3 4.258971 6.272973 2.418378 3.395009 5.340468 6.093108 4.324804 4.009715
## 4 3.983333 6.503333 2.145333 3.724733 5.192800 5.781933 4.280600 4.185222
## 5 3.982542 6.122935 1.942009 4.036441 5.664348 5.276085 4.193115 4.397148
## 6 4.378187 5.115919 2.175338 4.137515 5.404914 4.974350 4.235859 4.410030
## IriPul IriRei IriRuf LanAur LanFul LanLeu LanVer LopGri
## 1 2.556881 1.6641640 2.398438 1.711784 2.155905 6.806644 1.875807 4.396759
## 2 2.662247 1.6932523 2.305449 2.010827 2.309213 8.114988 1.914119 4.427661
## 3 2.329413 1.1432036 2.183699 1.496450 2.127450 8.692264 2.270701 4.317014
## 4 2.325067 1.1347333 2.418600 1.369267 1.809733 9.683662 2.384200 3.895667
## 5 2.009063 1.2780593 3.164674 1.538796 1.728348 10.244642 2.234641 3.910170
## 6 1.805400 0.9912414 3.542218 1.353515 1.455249 9.401372 2.153353 4.449196
## LopPus LoxAno LoxNoc LoxPor LoxVio MelMel MelNig MelRic
## 1 7.488234 3.910953 3.112941 0.7241694 3.978859 6.290421 3.317616 2.663009
## 2 7.728872 3.611694 3.302369 0.5207532 3.494323 6.346495 3.616895 2.822505
## 3 8.120056 3.903131 3.412677 0.4182378 2.825022 5.101267 3.853908 2.791955
## 4 8.777200 3.814917 3.506333 0.8708667 2.716667 4.902100 3.827133 1.614333
## 5 9.125893 3.826162 3.245481 0.7432741 2.436115 4.756217 4.466144 2.336722
## 6 9.077670 4.258948 2.971335 0.7415405 2.341905 4.632594 4.545243 2.238441
## MelXan NemPil NeoFas NepOne NesAcu OrcAbe OreFra OryAng
## 1 3.227230 4.226787 5.055903 2.750476 4.476550 4.686505 4.003340 4.710562
## 2 3.053802 4.427649 5.578859 2.984126 4.927500 5.367763 4.011671 4.578402
## 3 3.187752 5.991195 5.745822 3.277987 5.451937 6.685174 3.800905 4.211364
## 4 3.999750 6.571733 5.922933 4.008067 5.517333 5.427667 4.056417 3.453000
## 5 4.214352 6.875541 6.054052 3.772852 5.130167 5.959997 4.003097 3.298693
## 6 3.863957 7.280532 5.968306 3.802959 5.445279 5.704658 3.699477 3.357095
## OryCra OryFun OryMax OryNut ParCap ParCor ParDom ParGul
## 1 7.429500 3.171680 2.478050 3.547213 3.131365 2.192040 3.488714 5.262649
## 2 7.114906 1.696083 2.238608 3.181813 3.250432 1.692593 3.211316 5.341650
## 3 6.791849 1.783286 2.637338 3.065804 2.936459 1.480622 2.503977 5.598689
## 4 6.554433 2.014250 2.274000 3.493822 2.827667 1.769867 2.564167 5.632500
## 5 6.204209 1.687894 2.008000 2.797590 2.573741 2.164226 3.877600 5.221044
## 6 6.324417 2.340734 2.263856 3.435716 1.892505 2.509106 4.650066 5.927428
## ParHum ParNig PhrAla PhrAtr PhrCar PhrDor PhrEry PhrFru
## 1 1.291459 2.886041 9.293345 3.329982 4.506378 5.958768 5.393117 4.035439
## 2 1.470869 3.250788 9.722730 3.493220 4.852748 6.540127 5.173101 3.700482
## 3 2.591761 3.799027 9.971624 3.316380 5.439523 7.434704 5.215481 3.711192
## 4 3.333333 3.780667 10.363583 3.523467 5.290667 6.914625 5.228400 4.010433
## 5 1.288852 4.100083 11.004463 3.808385 6.077074 7.631400 5.338504 4.176483
## 6 2.063306 3.766527 11.143282 4.008009 5.946428 7.364779 5.950004 4.193376
## PhrPat PhrPle PhrPun PhrUni PieCin PinIno PipMel PlaCra
## 1 5.918407 4.492483 4.337694 4.419768 4.849804 7.595793 44.03170 4.935750
## 2 4.687977 4.608919 4.675339 4.146243 5.063029 7.981541 45.27175 4.996324
## 3 3.553551 4.704658 4.880333 4.382067 5.183631 6.137357 46.54217 5.160751
## 4 3.805467 4.587800 4.570267 4.332867 5.437083 5.480333 46.97613 5.289200
## 5 4.819904 4.487626 4.711126 3.984511 5.364856 4.238993 47.44355 5.450237
## 6 5.570063 4.528342 4.649456 4.118596 5.041036 3.369950 47.98022 5.174789
## PooAlt PooBol PooCab PooCae PooCin PooEry PooHis PooHyp
## 1 2.659149 5.074505 3.147468 1.979581 3.787056 7.712095 2.747193 4.714058
## 2 3.038234 4.900971 3.120582 2.135532 4.244384 7.609276 2.777488 4.660856
## 3 3.145950 4.980750 3.085348 2.241360 4.363910 7.682692 2.762346 4.884816
## 4 3.486833 5.398733 3.029067 2.709000 4.223600 7.945467 2.773867 5.114267
## 5 3.895370 5.529215 3.251222 3.093833 4.675441 8.458844 3.236433 4.961119
## 6 3.539613 5.523874 3.204924 2.964775 4.821710 8.432375 3.034631 4.787085
## PooMel PooNig PooOrn PooTho PooTor PooWhi PorCae PyrRuf
## 1 11.26016 5.275041 5.698532 4.047944 1.994281 6.175694 8.884247 2.985196
## 2 11.80131 5.415951 5.689101 4.329250 2.436692 6.643366 9.345823 2.704384
## 3 11.56038 5.513935 4.807838 4.291748 3.602207 6.608038 9.545780 3.434103
## 4 11.83930 4.856200 3.447583 4.118200 3.506933 6.874067 9.865067 3.132467
## 5 10.55111 3.877707 4.733398 4.123044 3.973230 6.915759 9.460619 4.185693
## 6 9.37007 4.818914 4.960957 4.444984 3.475429 6.550272 9.483804 3.727236
## RamBre RamCar RamCos RamDim RamFla RamMel RamNig
## 1 3.8178378 0.6575369 0.00000000 0.0000000 0.0000000 2.373670 0.8585081
## 2 1.7594474 1.2629982 0.00000000 0.4934586 0.0000000 2.212339 0.7717676
## 3 2.2941652 1.8017676 0.03942703 0.9421703 0.1915694 1.945461 1.1554324
## 4 0.7708889 2.4197333 0.28406667 1.0868333 0.0000000 1.657600 1.4081333
## 5 0.5567407 2.7865519 0.34956667 0.8162907 0.2692222 1.029870 1.4319370
## 6 2.3220991 2.3640991 0.31947207 0.3772622 0.0000000 1.078897 1.7661766
## RamPas RamSan RhoCru SalAlb SalAte SalAto SalAtp SalAur
## 1 2.092480 1.107279 5.267655 3.915486 0.00000000 4.418674 3.056477 4.420184
## 2 1.974804 1.321750 5.313646 3.759941 0.00000000 4.314685 2.939376 4.057112
## 3 1.886745 1.813571 5.048610 4.188431 0.00000000 4.326058 2.848426 4.018236
## 4 1.575917 2.131267 4.617833 3.978600 0.00000000 3.803733 2.996133 3.447733
## 5 1.236782 2.613637 4.703222 3.953630 0.04198148 3.729526 3.649130 3.590822
## 6 1.318827 2.669377 4.916980 4.746939 0.53640541 3.898820 3.974416 4.002832
## SalCin SalCoe SalFul SalGro SalMal SalMam SalMul SalNig
## 1 3.168025 4.004351 1.0794342 5.033497 5.002832 5.334402 15.74235 3.370532
## 2 2.839908 3.822681 1.0258180 4.819243 4.476708 5.543121 16.62963 3.105339
## 3 3.011948 3.490829 0.5417351 4.965391 4.330586 5.199097 15.40084 3.240399
## 4 3.607250 3.170400 2.0172667 4.799333 4.866000 4.940267 15.30760 3.271667
## 5 3.169894 3.324511 2.5264296 4.767463 6.317667 5.091256 14.95960 3.644216
## 6 3.618851 3.915580 2.5987027 4.726150 6.961162 5.189872 14.91429 4.024294
## SalOre SalRuf SalSim SalStr SchMel SchRuf SerAlb SicAur
## 1 5.463356 5.274351 4.474921 6.131899 4.324312 4.449703 14.00662 8.577032
## 2 7.085903 5.478189 4.156344 6.256351 4.736539 5.146381 14.34296 8.028575
## 3 6.254764 5.436554 4.238843 6.078312 4.933793 5.448093 14.93730 7.989236
## 4 5.441367 4.852167 4.596400 6.137133 4.674733 5.316222 15.42970 6.809067
## 5 5.154061 5.456111 4.878281 6.368189 5.512852 5.474654 16.32746 5.878826
## 6 4.032126 5.520279 5.259483 6.556189 5.823301 5.291928 16.89555 6.798753
## SicCit SicCol SicFla SicLeb SicLua SicLuc SicLuo SicLuv
## 1 11.25139 7.205077 2.830766 7.648045 4.676791 0.1021811 5.362960 5.264173
## 2 12.29999 7.287677 3.230378 7.145766 5.741486 0.2458739 5.011879 5.571210
## 3 12.99739 7.838106 4.078973 8.224568 6.527688 0.6680757 4.439155 4.883621
## 4 13.91120 6.471133 4.317333 7.508000 6.496067 1.2284667 4.784400 3.816833
## 5 14.09174 7.217593 4.263859 7.562741 5.551548 1.5280704 5.216737 3.680078
## 6 14.09012 6.273910 4.061829 7.912486 5.579627 1.6636676 5.659854 3.144914
## SicOli SicRai SicTac SicUro SpoAlb SpoAme SpoArd SpoBoe
## 1 5.204306 6.755234 4.564286 3.310614 4.268679 2.761505 3.879212 0.2185383
## 2 5.255086 6.567703 4.546519 3.668602 4.282187 3.020553 4.298396 1.2393041
## 3 5.441187 5.799820 4.824760 4.187177 4.633883 3.492620 4.489554 1.7047027
## 4 5.620867 6.934000 4.833533 4.504667 4.380400 3.838067 3.044667 1.8463333
## 5 5.716030 8.905722 5.345681 4.650185 4.569707 2.745678 2.185917 1.6927407
## 6 6.155941 8.194748 5.210966 4.254263 4.497658 2.635780 1.390338 1.2230000
## SpoBoo SpoCae SpoCas SpoCin SpoCol SpoCor SpoFro SpoHyc
## 1 4.993282 1.9897946 4.061829 6.635216 4.723308 4.821036 4.289338 10.106671
## 2 4.729613 1.4693568 3.707445 7.844453 4.328332 5.098705 3.536815 9.856275
## 3 4.072123 1.2005856 3.342933 8.342039 4.142750 7.214793 2.711018 10.156014
## 4 4.578778 0.9673333 3.272533 8.349444 4.217000 6.789800 3.588333 10.500500
## 5 5.105025 1.6696407 3.542422 8.371451 3.992930 7.721230 2.995056 10.530398
## 6 5.205201 1.9081766 3.553863 8.365844 4.160517 6.525811 3.909784 10.799360
## SpoHyx SpoIna SpoInt SpoLeu SpoLin SpoLuc SpoMin SpoNic
## 1 9.887982 11.759180 4.064547 13.144916 5.005912 4.448468 4.661175 5.567274
## 2 9.209042 10.563775 4.201367 12.932435 4.390182 4.509087 5.045897 5.589960
## 3 9.126309 9.028432 3.419466 12.430396 3.883852 4.325381 5.108537 5.457667
## 4 9.069533 9.184333 3.452333 10.863333 4.033333 4.926000 5.080867 5.027067
## 5 7.520901 9.913889 3.455208 9.510302 3.316396 4.462105 5.410693 5.467456
## 6 8.636017 11.347477 3.394881 9.717895 3.347627 4.186982 5.259308 6.098886
## SpoPer SpoPlu SpoRuf SpoSch SpoSim SpoTor SteDia TacCor
## 1 3.483946 5.242597 4.825178 5.069198 5.544849 3.914820 41.60987 4.836622
## 2 3.227615 5.732869 4.247578 5.417510 5.466721 4.728831 41.75434 5.258086
## 3 3.497984 6.309418 2.992063 5.246964 5.149845 4.435311 42.64125 5.863616
## 4 3.149917 6.454133 2.951133 5.253867 5.456667 3.899500 44.07067 6.663067
## 5 3.446389 6.809431 3.171059 5.042793 5.492134 3.539630 45.07088 6.291722
## 6 3.981626 6.210236 3.725485 4.817861 5.919448 2.797162 46.70900 6.223065
## TacCri TacDel TacLuc TacPho TacRuf TacRuv TacSur TanArg
## 1 1.640398 8.762058 3.119652 0.7513315 8.583155 2.088650 1.789730 2.427432
## 2 1.815859 8.307054 3.777277 1.2612685 7.230387 1.740654 1.856085 2.181986
## 3 1.351065 7.036616 4.345204 1.2830036 4.991850 2.443249 1.607371 3.051023
## 4 1.128133 7.506867 3.852400 0.8894667 3.590933 3.159267 1.926733 2.819000
## 5 1.416974 7.649233 3.349419 0.5677444 2.922400 3.260626 1.982904 3.043611
## 6 1.414360 7.682982 3.427677 0.4146937 3.382674 3.840914 1.735825 3.130450
## TanArt TanCal TanCay TanChi TanChr TanCuc TanCya TanCyc
## 1 1.467633 3.939865 2.766577 1.546768 1.677360 5.353996 12.85256 24.99214
## 2 1.070534 4.167627 3.157926 2.042294 1.524777 4.903018 12.69680 25.51809
## 3 1.602883 4.549110 3.320213 2.138168 1.414596 5.572593 13.72460 26.39848
## 4 1.291167 4.865400 3.444867 2.011733 1.588752 6.040467 14.50273 26.02467
## 5 1.602079 5.802404 3.449726 2.003015 1.498061 5.336267 14.47962 26.77763
## 6 1.876052 6.399480 3.999852 2.592816 1.696240 6.013441 15.54088 28.33923
## TanCyo TanCyp TanCyv TanDes TanDow TanFas TanFlo TanFuc
## 1 4.685744 4.894128 1.53088288 3.585494 8.093189 7.288550 2.364691 4.202508
## 2 3.982295 5.122881 0.39092793 3.767771 7.414797 8.486483 2.275387 4.517622
## 3 4.159786 5.076450 0.01966066 3.400597 8.021703 9.367165 2.907186 4.599910
## 4 3.679571 5.279167 0.00000000 2.275361 7.958667 9.500778 3.105667 4.445889
## 5 3.799271 5.238824 0.00000000 2.683393 9.429824 8.347049 2.306123 5.335204
## 6 4.048019 5.988795 0.00000000 2.766825 11.963306 7.516703 2.475210 5.514745
## TanGut TanGyr TanHei TanIct TanIno TanJoh TanLab TanLar
## 1 4.829768 3.375108 4.754764 1.661829 4.501844 5.236825 8.742186 7.501444
## 2 4.057305 3.741441 4.464932 1.760623 4.582871 5.657435 10.544344 7.244984
## 3 3.251351 3.843132 3.552948 2.138634 4.196477 6.356735 9.801495 7.059277
## 4 3.480000 4.533933 3.714467 2.309267 4.342778 5.917067 10.277800 7.178000
## 5 3.814641 4.747278 1.882696 2.231774 4.564636 5.688391 9.468952 7.093171
## 6 3.383431 4.935061 1.599903 2.125413 4.895985 6.643610 9.604049 7.036554
## TanLav TanMex TanMey TanNic TanNiv TanPal TanPar TanPer
## 1 2.290441 2.827105 2.478378 21.84156 8.460319 11.34139 4.762243 5.691432
## 2 2.309212 3.017802 2.591824 22.41399 8.460668 11.56592 4.895393 6.298126
## 3 2.189450 2.984490 3.343730 23.14435 8.719510 12.82801 5.003978 6.256604
## 4 1.683000 3.178200 3.253833 23.55653 8.946267 13.70775 5.038467 6.772000
## 5 1.570694 3.372519 3.030491 24.11456 8.996733 13.53459 5.232363 6.668056
## 6 2.021892 3.450267 3.097725 25.70736 9.022232 14.43345 5.463470 5.572568
## TanPre TanPun TanRue TanRuf TanRuu TanSch TanSel TanVar
## 1 9.944220 6.008441 5.002063 3.424358 6.909399 1.863059 10.139715 3.795550
## 2 9.816647 6.111986 5.386405 2.702604 6.704295 1.906650 9.674646 3.480396
## 3 10.643818 6.772058 4.987523 2.554797 6.862164 1.649207 9.437048 3.359012
## 4 10.704467 7.040267 4.996000 2.483500 6.317833 1.751800 10.165000 3.968000
## 5 10.522515 7.286067 5.009926 2.402514 6.948602 1.717885 10.715284 4.519877
## 6 10.092454 7.242719 5.891712 2.773025 6.433117 1.826332 11.083459 5.116051
## TanVas TanVel TanVir TanVit TanXac TanXag TerVir ThlFul
## 1 2.090986 3.547241 5.851671 0.0000000 3.319324 5.198502 22.83156 6.217622
## 2 1.932000 3.829730 5.830550 0.0000000 3.283024 5.814634 23.48568 5.804049
## 3 1.741870 3.966005 5.949300 0.3065514 3.626273 6.880012 23.79442 6.021079
## 4 2.395533 4.065520 6.169000 0.3940000 3.822667 7.266778 25.07140 6.482333
## 5 3.007759 4.212550 6.293051 0.7635556 3.874549 7.289660 25.85069 6.294207
## 6 3.461418 4.419631 6.829633 0.8666216 4.072607 7.429174 26.23603 6.344178
## ThlIno ThlOrn ThlPec ThlRuf ThlSor ThrAbb ThrBon ThrCya
## 1 4.137212 4.018599 4.602153 4.298849 2.225132 16.54398 17.06939 23.94336
## 2 4.004766 4.506939 5.335471 4.226583 2.187126 16.95915 17.37653 23.97398
## 3 4.393178 5.093394 5.356084 4.089910 2.317160 17.08833 18.18349 24.66325
## 4 3.961167 5.019833 5.259889 3.864750 2.380533 17.96227 19.00967 26.03373
## 5 3.180060 4.355625 5.754611 3.825542 2.680167 18.30083 19.33483 26.94218
## 6 3.277047 4.366608 4.740420 3.410989 3.460623 18.89690 19.63418 27.52753
## ThrCyp ThrEpi ThrGla ThrOrn ThrPal ThrSay TiaBic TiaCan
## 1 11.77104 15.70287 11.164295 26.62071 11.45043 7.734886 0.8755788 5.133854
## 2 12.07851 15.92114 12.126169 26.56315 11.98017 7.490996 1.3769189 5.877162
## 3 11.92041 16.61791 11.888736 27.70014 12.96002 8.468750 2.1457500 6.069865
## 4 10.88247 17.18853 10.780833 29.30667 13.44800 9.263733 3.4594167 5.848733
## 5 11.06614 17.83998 9.608986 28.68717 13.98219 9.005744 3.1475093 5.876067
## 6 10.79618 17.93138 8.779264 28.67975 14.47977 9.377237 2.8914527 5.991153
## TiaFul TiaObs TiaOli TriMel UroSto VolJac WetSte XenCon
## 1 1.579715 4.562836 2.961937 4.663139 4.309291 4.797721 2.685617 6.996423
## 2 2.509440 4.996505 2.660923 4.703391 4.823276 4.793861 2.610338 7.475644
## 3 3.427998 6.191164 3.325404 5.140431 4.588535 5.222081 3.710387 6.674491
## 4 4.815200 5.576417 3.770933 5.100200 4.847111 5.361533 4.143833 6.251667
## 5 4.468637 6.331551 3.296244 4.975270 5.256043 5.136441 3.369602 5.716065
## 6 4.154216 6.832324 3.227438 4.348146 5.330523 4.315382 3.279829 5.451829
## XenPar
## 1 11.79418
## 2 12.58418
## 3 12.37302
## 4 12.37817
## 5 13.14826
## 6 13.42591Step 1 Extracting tanager tree
get the tanagers only
treepath <- "/Users/michaelalfaro/Dropbox/color_evolution_scratch/tanager_comp_analysis/allison_data/MCC_Tree_SpNames.nex"
tt <- read.nexus(treepath)
# Extract the tip labels (species names) from the tree
tan_species <- grep("^Tan", tt$tip.label, value = TRUE)
# Find the MRCA node of the species that start with "Tan"
mrca_node <- getMRCA(tt, tan_species)
# Extract the subtree from the MRCA node
tanagerTree <- extract.clade(tt, mrca_node)
# Plot the tanagerTree
p <- ggtree(tanagerTree) +
geom_tiplab(size = 2.5, align = TRUE) +
theme_tree2()
# Display the plot
print(p)
### Step 2 Create a Data Frame Matching the Tanager Tree
```r
# Extract the species names (tip labels) from the tanagerTree
tanager_species <- tanagerTree$tip.label
# Filter crown_data to include only species that are present in tanagerTree
tanData_filtered <- crown_data %>%
select(wl, any_of(tanager_species))
# Display the filtered tanData to verify
head(tanData_filtered)## wl TanCyo TanLab TanRue TanChi TanVel TanCal TanMex TanIno
## 1 300 4.685744 8.742186 5.002063 1.546768 3.547241 3.939865 2.827105 4.501844
## 2 301 3.982295 10.544344 5.386405 2.042294 3.829730 4.167627 3.017802 4.582871
## 3 302 4.159786 9.801495 4.987523 2.138168 3.966005 4.549110 2.984490 4.196477
## 4 303 3.679571 10.277800 4.996000 2.011733 4.065520 4.865400 3.178200 4.342778
## 5 304 3.799271 9.468952 5.009926 2.003015 4.212550 5.802404 3.372519 4.564636
## 6 305 4.048019 9.604049 5.891712 2.592816 4.419631 6.399480 3.450267 4.895985
## TanSel TanFas TanCyc TanCyv TanDes TanChr TanXac TanSch
## 1 10.139715 7.288550 24.99214 1.53088288 3.585494 1.677360 3.319324 1.863059
## 2 9.674646 8.486483 25.51809 0.39092793 3.767771 1.524777 3.283024 1.906650
## 3 9.437048 9.367165 26.39848 0.01966066 3.400597 1.414596 3.626273 1.649207
## 4 10.165000 9.500778 26.02467 0.00000000 2.275361 1.588752 3.822667 1.751800
## 5 10.715284 8.347049 26.77763 0.00000000 2.683393 1.498061 3.874549 1.717885
## 6 11.083459 7.516703 28.33923 0.00000000 2.766825 1.696240 4.072607 1.826332
## TanJoh TanArt TanIct TanFlo TanPar TanLav TanGyr TanVas
## 1 5.236825 1.467633 1.661829 2.364691 4.762243 2.290441 3.375108 2.090986
## 2 5.657435 1.070534 1.760623 2.275387 4.895393 2.309212 3.741441 1.932000
## 3 6.356735 1.602883 2.138634 2.907186 5.003978 2.189450 3.843132 1.741870
## 4 5.917067 1.291167 2.309267 3.105667 5.038467 1.683000 4.533933 2.395533
## 5 5.688391 1.602079 2.231774 2.306123 5.232363 1.570694 4.747278 3.007759
## 6 6.643610 1.876052 2.125413 2.475210 5.463470 2.021892 4.935061 3.461418
## TanNiv TanFuc TanDow TanRuu TanPun TanXag TanGut TanVar
## 1 8.460319 4.202508 8.093189 6.909399 6.008441 5.198502 4.829768 3.795550
## 2 8.460668 4.517622 7.414797 6.704295 6.111986 5.814634 4.057305 3.480396
## 3 8.719510 4.599910 8.021703 6.862164 6.772058 6.880012 3.251351 3.359012
## 4 8.946267 4.445889 7.958667 6.317833 7.040267 7.266778 3.480000 3.968000
## 5 8.996733 5.335204 9.429824 6.948602 7.286067 7.289660 3.814641 4.519877
## 6 9.022232 5.514745 11.963306 6.433117 7.242719 7.429174 3.383431 5.116051
## TanPal TanVit TanCuc TanCay TanPre TanNic TanCya TanLar
## 1 11.34139 0.0000000 5.353996 2.766577 9.944220 21.84156 12.85256 7.501444
## 2 11.56592 0.0000000 4.903018 3.157926 9.816647 22.41399 12.69680 7.244984
## 3 12.82801 0.3065514 5.572593 3.320213 10.643818 23.14435 13.72460 7.059277
## 4 13.70775 0.3940000 6.040467 3.444867 10.704467 23.55653 14.50273 7.178000
## 5 13.53459 0.7635556 5.336267 3.449726 10.522515 24.11456 14.47962 7.093171
## 6 14.43345 0.8666216 6.013441 3.999852 10.092454 25.70736 15.54088 7.036554
## TanCyp TanVir TanHei TanArg ThrCyp ThrAbb ThrPal ThrOrn
## 1 4.894128 5.851671 4.754764 2.427432 11.77104 16.54398 11.45043 26.62071
## 2 5.122881 5.830550 4.464932 2.181986 12.07851 16.95915 11.98017 26.56315
## 3 5.076450 5.949300 3.552948 3.051023 11.92041 17.08833 12.96002 27.70014
## 4 5.279167 6.169000 3.714467 2.819000 10.88247 17.96227 13.44800 29.30667
## 5 5.238824 6.293051 1.882696 3.043611 11.06614 18.30083 13.98219 28.68717
## 6 5.988795 6.829633 1.599903 3.130450 10.79618 18.89690 14.47977 28.67975
## ThrSay ThrEpi TanRuf
## 1 7.734886 15.70287 3.424358
## 2 7.490996 15.92114 2.702604
## 3 8.468750 16.61791 2.554797
## 4 9.263733 17.18853 2.483500
## 5 9.005744 17.83998 2.402514
## 6 9.377237 17.93138 2.773025Step 3: Fitting Natural and Cubic Splines to the Tanager Data
### Step 3: Fitting Natural and Cubic Splines to the Tanager Data
# Function to fit a natural spline and a cubic spline, returning the models, coefficients, goodness-of-fit metrics, and RGB colors
fit_splines <- function(wl, reflectance) {
# Fit natural cubic spline
natural_spline_fit <- lm(reflectance ~ ns(wl, df = 5))
# Fit cubic spline
cubic_spline_fit <- lm(reflectance ~ poly(wl, 3, raw = TRUE))
# Convert fitted splines to rspec and then to RGB hex codes
natural_spline_rspec <- as.rspec(data.frame(wl = wl, reflectance = predict(natural_spline_fit)), lim = c(300, 700))
natural_rgb <- spec2rgb(natural_spline_rspec)
cubic_spline_rspec <- as.rspec(data.frame(wl = wl, reflectance = predict(cubic_spline_fit)), lim = c(300, 700))
cubic_rgb <- spec2rgb(cubic_spline_rspec)
# Extract coefficients
natural_coef <- round(coef(natural_spline_fit), 2)
cubic_coef <- round(coef(cubic_spline_fit), 2)
# Calculate goodness-of-fit metrics
natural_adj_r2 <- summary(natural_spline_fit)$adj.r.squared
cubic_adj_r2 <- summary(cubic_spline_fit)$adj.r.squared
natural_aic <- AIC(natural_spline_fit)
cubic_aic <- AIC(cubic_spline_fit)
return(list(
natural_spline_fit = natural_spline_fit,
cubic_spline_fit = cubic_spline_fit,
natural_rgb = natural_rgb,
cubic_rgb = cubic_rgb,
natural_coef = natural_coef,
cubic_coef = cubic_coef,
natural_adj_r2 = natural_adj_r2,
cubic_adj_r2 = cubic_adj_r2,
natural_aic = natural_aic,
cubic_aic = cubic_aic
))
}
# Apply the spline fitting function to each species
spline_results <- lapply(tanData_filtered[-1], function(reflectance) {
fit_splines(tanData_filtered$wl, reflectance)
})## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1## The spectral data contain 82 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 82 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : 1 negative quantum-catch value(s) returned, which may be unreliable.
## Consider whether the illuminant is properly scaled, and the appropriate form of
## quantum catch is being calculated.## wavelengths found in column 1## The spectral data contain 89 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 89 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1## The spectral data contain 15 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 15 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1## The spectral data contain 14 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 14 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1## The spectral data contain 37 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 37 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## The spectral data contain 8 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 8 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1## The spectral data contain 21 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 21 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## The spectral data contain 25 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 25 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## The spectral data contain 59 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 59 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1## The spectral data contain 77 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 77 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## The spectral data contain 27 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 27 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## The spectral data contain 48 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 48 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## The spectral data contain 63 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 63 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : 1 negative quantum-catch value(s) returned, which may be unreliable.
## Consider whether the illuminant is properly scaled, and the appropriate form of
## quantum catch is being calculated.## wavelengths found in column 1## The spectral data contain 78 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 78 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1## The spectral data contain 35 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 35 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1## The spectral data contain 22 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 22 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1## The spectral data contain 1 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 1 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1## The spectral data contain 30 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 30 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## The spectral data contain 3 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 3 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1## The spectral data contain 5 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 5 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1Step 3b–making the species plots of natural and cubic spline
plot_splines_with_metrics <- function(species_name, wl, reflectance, spline_results, rgb_colors) {
result <- spline_results[[species_name]]
par(mfrow = c(1, 2))
# Natural Spline
plot(wl, reflectance, type = "p", pch = 16, cex = 0.5, # Adjusted size of points
col = adjustcolor("black", alpha.f = 0.5), # Adjusted transparency of points
main = paste("Natural Spline -", species_name),
xlab = "Wavelength (nm)", ylab = "Reflectance")
lines(wl, predict(result$natural_spline_fit), col = rgb_colors[[species_name]], lwd = 2)
mtext(paste("Adj R2:", round(result$natural_adj_r2, 3), "AIC:", round(result$natural_aic, 2)),
side = 3, line = 0.5, cex = 0.8, col = "blue")
# Cubic Spline
plot(wl, reflectance, type = "p", pch = 16, cex = 0.5, # Adjusted size of points
col = adjustcolor("black", alpha.f = 0.5), # Adjusted transparency of points
main = paste("Cubic Spline -", species_name),
xlab = "Wavelength (nm)", ylab = "Reflectance")
lines(wl, predict(result$cubic_spline_fit), col = rgb_colors[[species_name]], lwd = 2)
mtext(paste("Adj R2:", round(result$cubic_adj_r2, 3), "AIC:", round(result$cubic_aic, 2)),
side = 3, line = 0.5, cex = 0.8, col = "blue")
}
for (species_name in names(spline_results)) {
plot_splines_with_metrics(species_name, tanData_filtered$wl, tanData_filtered[[species_name]], spline_results, rgb_colors)
}
pdf("natural_vs_cubic_splines_with_metrics.pdf", width = 12, height = 8)
for (species_name in names(spline_results)) {
plot_splines_with_metrics(species_name, tanData_filtered$wl, tanData_filtered[[species_name]], spline_results, rgb_colors)
}
dev.off()## quartz_off_screen
## 2Step 4: Plotting Spline Fits and Original Data in the Markdown Document
### Step 4: Calculate Average Splines and Standard Deviation
# Calculate the maximum reflectance across all splines
max_reflectance <- max(sapply(spline_results, function(result) {
max(predict(result$natural_spline_fit), predict(result$cubic_spline_fit))
}))
# Initialize vectors to store average splines and standard deviations
natural_spline_avg <- numeric(length(tanData_filtered$wl))
cubic_spline_avg <- numeric(length(tanData_filtered$wl))
natural_spline_sd <- numeric(length(tanData_filtered$wl))
cubic_spline_sd <- numeric(length(tanData_filtered$wl))
# Calculate averages
for (result in spline_results) {
natural_pred <- predict(result$natural_spline_fit)
cubic_pred <- predict(result$cubic_spline_fit)
natural_spline_avg <- natural_spline_avg + natural_pred
cubic_spline_avg <- cubic_spline_avg + cubic_pred
}
# Convert sums to averages
natural_spline_avg <- natural_spline_avg / length(spline_results)
cubic_spline_avg <- cubic_spline_avg / length(spline_results)
# Calculate standard deviations
for (result in spline_results) {
natural_pred <- predict(result$natural_spline_fit)
cubic_pred <- predict(result$cubic_spline_fit)
natural_spline_sd <- natural_spline_sd + (natural_pred - natural_spline_avg)^2
cubic_spline_sd <- cubic_spline_sd + (cubic_pred - cubic_spline_avg)^2
}
natural_spline_sd <- sqrt(natural_spline_sd / (length(spline_results) - 1))
cubic_spline_sd <- sqrt(cubic_spline_sd / (length(spline_results) - 1))
### Correct for negative values in the average splines
# Convert the average natural and cubic splines to rspec objects
natural_spline_rspec <- as.rspec(data.frame(wl = tanData_filtered$wl, reflectance = natural_spline_avg), lim = c(300, 700))## wavelengths found in column 1cubic_spline_rspec <- as.rspec(data.frame(wl = tanData_filtered$wl, reflectance = cubic_spline_avg), lim = c(300, 700))## wavelengths found in column 1# Correct any negative values in the rspec objects
natural_spline_rspec <- procspec(natural_spline_rspec, fixneg = "zero")## processing options applied:
## Negative value correction: converted negative values to zerocubic_spline_rspec <- procspec(cubic_spline_rspec, fixneg = "zero")## processing options applied:
## Negative value correction: converted negative values to zero# Convert corrected rspec objects to RGB hex codes
natural_spline_rgb <- spec2rgb(natural_spline_rspec)
cubic_spline_rgb <- spec2rgb(cubic_spline_rspec)
### Step 4a: Plot Average Splines with Confidence Intervals
# Create a plot for the average natural spline with confidence intervals
plot(tanData_filtered$wl, natural_spline_avg, type = "l", col = natural_spline_rgb, lwd = 4,
main = "Average Natural Spline with Variation Hotspots",
xlab = "Wavelength (nm)", ylab = "Reflectance", ylim = c(0, max_reflectance))
polygon(c(tanData_filtered$wl, rev(tanData_filtered$wl)),
c(natural_spline_avg + natural_spline_sd, rev(natural_spline_avg - natural_spline_sd)),
col = adjustcolor(natural_spline_rgb, alpha.f = 0.2), border = NA)
lines(tanData_filtered$wl, natural_spline_avg, col = natural_spline_rgb, lwd = 4)
# Create a plot for the average cubic spline with confidence intervals
plot(tanData_filtered$wl, cubic_spline_avg, type = "l", col = cubic_spline_rgb, lwd = 4,
main = "Average Cubic Spline with Variation Hotspots",
xlab = "Wavelength (nm)", ylab = "Reflectance", ylim = c(0, max_reflectance))
polygon(c(tanData_filtered$wl, rev(tanData_filtered$wl)),
c(cubic_spline_avg + cubic_spline_sd, rev(cubic_spline_avg - cubic_spline_sd)),
col = adjustcolor(cubic_spline_rgb, alpha.f = 0.2), border = NA)
lines(tanData_filtered$wl, cubic_spline_avg, col = cubic_spline_rgb, lwd = 4)
# Save the hotspot plots to a PDF
pdf("spline_hotspots.pdf", width = 12, height = 8)
par(mfrow = c(1, 2)) # Plot both splines side by side
# Average Natural Spline with Hotspots in PDF
plot(tanData_filtered$wl, natural_spline_avg, type = "l", col = natural_spline_rgb, lwd = 4,
main = "Average Natural Spline with Variation Hotspots",
xlab = "Wavelength (nm)", ylab = "Reflectance", ylim = c(0, max_reflectance))
polygon(c(tanData_filtered$wl, rev(tanData_filtered$wl)),
c(natural_spline_avg + natural_spline_sd, rev(natural_spline_avg - natural_spline_sd)),
col = adjustcolor(natural_spline_rgb, alpha.f = 0.2), border = NA)
lines(tanData_filtered$wl, natural_spline_avg, col = natural_spline_rgb, lwd = 4)
# Average Cubic Spline with Hotspots in PDF
plot(tanData_filtered$wl, cubic_spline_avg, type = "l", col = cubic_spline_rgb, lwd = 4,
main = "Average Cubic Spline with Variation Hotspots",
xlab = "Wavelength (nm)", ylab = "Reflectance", ylim = c(0, max_reflectance))
polygon(c(tanData_filtered$wl, rev(tanData_filtered$wl)),
c(cubic_spline_avg + cubic_spline_sd, rev(cubic_spline_avg - cubic_spline_sd)),
col = adjustcolor(cubic_spline_rgb, alpha.f = 0.2), border = NA)
lines(tanData_filtered$wl, cubic_spline_avg, col = cubic_spline_rgb, lwd = 4)
dev.off()## quartz_off_screen
## 2### Step 4b: Plot Natural Spline Mean and Cubic Spline Mean with Species Overlaid
# Plot Natural Spline Mean with all species overlaid
par(mfrow = c(1, 2))
# Plot Natural Spline Mean with Species Overlaid
plot(tanData_filtered$wl, natural_spline_avg, type = "l", col = natural_spline_rgb, lwd = 4,
main = "Natural Spline Mean with Species Overlaid",
xlab = "Wavelength (nm)", ylab = "Reflectance", ylim = c(0, max_reflectance))
for (i in seq_along(spline_results)) {
species_rspec <- as.rspec(data.frame(wl = tanData_filtered$wl, reflectance = predict(spline_results[[i]]$natural_spline_fit)), lim = c(300, 700))
species_rspec <- procspec(species_rspec, fixneg = "zero") # Correct any negative values
species_rgb <- spec2rgb(species_rspec)
lines(tanData_filtered$wl, predict(spline_results[[i]]$natural_spline_fit), col = adjustcolor(species_rgb, alpha.f = 0.3))
}## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 82 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 37 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 25 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 27 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 63 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 30 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zerolines(tanData_filtered$wl, natural_spline_avg, col = natural_spline_rgb, lwd = 4)
# Plot Cubic Spline Mean with Species Overlaid
plot(tanData_filtered$wl, cubic_spline_avg, type = "l", col = cubic_spline_rgb, lwd = 4,
main = "Cubic Spline Mean with Species Overlaid",
xlab = "Wavelength (nm)", ylab = "Reflectance", ylim = c(0, max_reflectance))
for (i in seq_along(spline_results)) {
species_rspec <- as.rspec(data.frame(wl = tanData_filtered$wl, reflectance = predict(spline_results[[i]]$cubic_spline_fit)), lim = c(300, 700))
species_rspec <- procspec(species_rspec, fixneg = "zero") # Correct any negative values
species_rgb <- spec2rgb(species_rspec)
lines(tanData_filtered$wl, predict(spline_results[[i]]$cubic_spline_fit), col = adjustcolor(species_rgb, alpha.f = 0.3))
}## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 89 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 15 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 14 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 8 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 21 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 59 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 77 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 48 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 78 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 35 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 22 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 1 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 3 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 5 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zerolines(tanData_filtered$wl, cubic_spline_avg, col = cubic_spline_rgb, lwd = 4)
# Save the overlay plots to a PDF
pdf("spline_mean_with_species_overlay.pdf", width = 12, height = 8)
par(mfrow = c(1, 2))
# Plot Natural Spline Mean with Species Overlaid in PDF
plot(tanData_filtered$wl, natural_spline_avg, type = "l", col = natural_spline_rgb, lwd = 4,
main = "Natural Spline Mean with Species Overlaid",
xlab = "Wavelength (nm)", ylab = "Reflectance", ylim = c(0, max_reflectance))
for (i in seq_along(spline_results)) {
species_rspec <- as.rspec(data.frame(wl = tanData_filtered$wl, reflectance = predict(spline_results[[i]]$natural_spline_fit)), lim = c(300, 700))
species_rspec <- procspec(species_rspec, fixneg = "zero") # Correct any negative values
species_rgb <- spec2rgb(species_rspec)
lines(tanData_filtered$wl, predict(spline_results[[i]]$natural_spline_fit), col = adjustcolor(species_rgb, alpha.f = 0.3))
}## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 82 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 37 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 25 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 27 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 63 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 30 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zerolines(tanData_filtered$wl, natural_spline_avg, col = natural_spline_rgb, lwd = 4)
# Plot Cubic Spline Mean with Species Overlaid in PDF
plot(tanData_filtered$wl, cubic_spline_avg, type = "l", col = cubic_spline_rgb, lwd = 4,
main = "Cubic Spline Mean with Species Overlaid",
xlab = "Wavelength (nm)", ylab = "Reflectance", ylim = c(0, max_reflectance))
for (i in seq_along(spline_results)) {
species_rspec <- as.rspec(data.frame(wl = tanData_filtered$wl, reflectance = predict(spline_results[[i]]$cubic_spline_fit)), lim = c(300, 700))
species_rspec <- procspec(species_rspec, fixneg = "zero") # Correct any negative values
species_rgb <- spec2rgb(species_rspec)
lines(tanData_filtered$wl, predict(spline_results[[i]]$cubic_spline_fit), col = adjustcolor(species_rgb, alpha.f = 0.3))
}## wavelengths found in column 1
## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 89 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 15 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 14 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 8 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 21 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 59 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 77 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 48 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 78 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 35 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 22 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 1 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 3 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## The spectral data contain 5 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zero## wavelengths found in column 1## processing options applied:
## Negative value correction: converted negative values to zerolines(tanData_filtered$wl, cubic_spline_avg, col = cubic_spline_rgb, lwd = 4)
dev.off()## quartz_off_screen
## 2Step 4b Make a pair of plots of the natural and cubic spline where the spline is colored as a spec that is converted to its hsv value (make sure no negative values)
### Step 4b: Plot Natural Spline with Species Overlaid and Print to PDF
# Convert the natural spline to an rspec object and then to an RGB hex code
natural_spline_rspec <- as.rspec(data.frame(wl = tanData_filtered$wl, reflectance = natural_spline_avg), lim = c(300, 700))## wavelengths found in column 1natural_spline_rgb <- spec2rgb(natural_spline_rspec)
# Convert the cubic spline to an rspec object and then to an RGB hex code
cubic_spline_rspec <- as.rspec(data.frame(wl = tanData_filtered$wl, reflectance = cubic_spline_avg), lim = c(300, 700))## wavelengths found in column 1cubic_spline_rgb <- spec2rgb(cubic_spline_rspec)
# Plot the natural and cubic splines side by side in R Markdown
par(mfrow = c(1, 2)) # Arrange plots side by side
# Plot Natural Spline with Species Overlaid
plot(tanData_filtered$wl, natural_spline_avg, type = "l", col = natural_spline_rgb, lwd = 4, # Increased thickness
main = "Natural Spline with Species Overlaid",
xlab = "Wavelength (nm)", ylab = "Reflectance", ylim = c(0, max_reflectance))
# Loop through each species to plot their splines with slightly more visible transparency
for (species in names(spline_results)) {
species_spline_rspec <- as.rspec(data.frame(wl = tanData_filtered$wl, reflectance = predict(spline_results[[species]]$natural_spline_fit)), lim = c(300, 700))
species_rgb <- spec2rgb(species_spline_rspec)
# Adjust transparency to make species lines slightly more visible
species_rgb_alpha <- adjustcolor(species_rgb, alpha.f = 0.4) # Adjusted alpha for more visibility
lines(tanData_filtered$wl, predict(spline_results[[species]]$natural_spline_fit), col = species_rgb_alpha)
}## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1## The spectral data contain 82 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 82 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : 1 negative quantum-catch value(s) returned, which may be unreliable.
## Consider whether the illuminant is properly scaled, and the appropriate form of
## quantum catch is being calculated.## wavelengths found in column 1## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1## The spectral data contain 37 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 37 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1
## wavelengths found in column 1## The spectral data contain 25 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 25 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1
## wavelengths found in column 1## The spectral data contain 27 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 27 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## The spectral data contain 63 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 63 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : 1 negative quantum-catch value(s) returned, which may be unreliable.
## Consider whether the illuminant is properly scaled, and the appropriate form of
## quantum catch is being calculated.## wavelengths found in column 1## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1## The spectral data contain 30 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 30 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1# Re-plot the natural spline on top to ensure it stands out
lines(tanData_filtered$wl, natural_spline_avg, col = natural_spline_rgb, lwd = 4) # Increased thickness
# Plot Cubic Spline with Species Overlaid
plot(tanData_filtered$wl, cubic_spline_avg, type = "l", col = cubic_spline_rgb, lwd = 4, # Increased thickness
main = "Cubic Spline with Species Overlaid",
xlab = "Wavelength (nm)", ylab = "Reflectance", ylim = c(0, max_reflectance))
# Loop through each species to plot their splines with slightly more visible transparency
for (species in names(spline_results)) {
species_spline_rspec <- as.rspec(data.frame(wl = tanData_filtered$wl, reflectance = predict(spline_results[[species]]$cubic_spline_fit)), lim = c(300, 700))
species_rgb <- spec2rgb(species_spline_rspec)
# Adjust transparency to make species lines slightly more visible
species_rgb_alpha <- adjustcolor(species_rgb, alpha.f = 0.4) # Adjusted alpha for more visibility
lines(tanData_filtered$wl, predict(spline_results[[species]]$cubic_spline_fit), col = species_rgb_alpha)
}## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1## The spectral data contain 89 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 89 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1## The spectral data contain 15 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 15 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## The spectral data contain 14 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 14 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1
## wavelengths found in column 1## The spectral data contain 8 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 8 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1## The spectral data contain 21 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 21 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## The spectral data contain 59 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 59 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1## The spectral data contain 77 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 77 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## The spectral data contain 48 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 48 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## The spectral data contain 78 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 78 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## The spectral data contain 35 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 35 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1
## wavelengths found in column 1## The spectral data contain 22 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 22 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## The spectral data contain 1 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 1 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1## The spectral data contain 3 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 3 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1## The spectral data contain 5 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 5 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1# Re-plot the cubic spline on top to ensure it stands out
lines(tanData_filtered$wl, cubic_spline_avg, col = cubic_spline_rgb, lwd = 4) # Increased thickness
# Save the plots to a PDF as well
pdf("natural_and_cubic_splines_with_species_overlay.pdf", width = 12, height = 8)
par(mfrow = c(1, 2)) # Arrange plots side by side in the PDF
# Plot Natural Spline with Species Overlaid in PDF
plot(tanData_filtered$wl, natural_spline_avg, type = "l", col = natural_spline_rgb, lwd = 4,
main = "Natural Spline with Species Overlaid",
xlab = "Wavelength (nm)", ylab = "Reflectance", ylim = c(0, max_reflectance))
for (species in names(spline_results)) {
species_rgb_alpha <- adjustcolor(spec2rgb(as.rspec(data.frame(wl = tanData_filtered$wl, reflectance = predict(spline_results[[species]]$natural_spline_fit)), lim = c(300, 700))), alpha.f = 0.4)
lines(tanData_filtered$wl, predict(spline_results[[species]]$natural_spline_fit), col = species_rgb_alpha)
}## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1## The spectral data contain 82 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 82 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : 1 negative quantum-catch value(s) returned, which may be unreliable.
## Consider whether the illuminant is properly scaled, and the appropriate form of
## quantum catch is being calculated.## wavelengths found in column 1## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1## The spectral data contain 37 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 37 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1
## wavelengths found in column 1## The spectral data contain 25 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 25 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1
## wavelengths found in column 1## The spectral data contain 27 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 27 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## The spectral data contain 63 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 63 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.
## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : 1 negative quantum-catch value(s) returned, which may be unreliable.
## Consider whether the illuminant is properly scaled, and the appropriate form of
## quantum catch is being calculated.## wavelengths found in column 1## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1## The spectral data contain 30 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 30 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1lines(tanData_filtered$wl, natural_spline_avg, col = natural_spline_rgb, lwd = 4)
# Plot Cubic Spline with Species Overlaid in PDF
plot(tanData_filtered$wl, cubic_spline_avg, type = "l", col = cubic_spline_rgb, lwd = 4,
main = "Cubic Spline with Species Overlaid",
xlab = "Wavelength (nm)", ylab = "Reflectance", ylim = c(0, max_reflectance))
for (species in names(spline_results)) {
species_rgb_alpha <- adjustcolor(spec2rgb(as.rspec(data.frame(wl = tanData_filtered$wl, reflectance = predict(spline_results[[species]]$cubic_spline_fit)), lim = c(300, 700))), alpha.f = 0.4)
lines(tanData_filtered$wl, predict(spline_results[[species]]$cubic_spline_fit), col = species_rgb_alpha)
}## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1## The spectral data contain 89 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 89 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1## The spectral data contain 15 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 15 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## The spectral data contain 14 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 14 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1
## wavelengths found in column 1## The spectral data contain 8 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 8 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1## The spectral data contain 21 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 21 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## The spectral data contain 59 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 59 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1## The spectral data contain 77 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 77 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## The spectral data contain 48 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 48 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## The spectral data contain 78 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 78 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## The spectral data contain 35 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 35 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1
## wavelengths found in column 1## The spectral data contain 22 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 22 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## The spectral data contain 1 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 1 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1## The spectral data contain 3 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 3 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1## The spectral data contain 5 negative value(s),
## which may produce unexpected results if used in models.
## Consider using procspec() to correct them.## Warning in vismodel(rspecdata, visual = "cie10", illum = "D65", vonkries =
## TRUE, : The spectral data contain 5 negative value(s), which may produce
## unexpected results. Consider using procspec() to correct them.## wavelengths found in column 1## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1
## wavelengths found in column 1lines(tanData_filtered$wl, cubic_spline_avg, col = cubic_spline_rgb, lwd = 4)
dev.off()## quartz_off_screen
## 2Step 5 Print a Summary Table of Spline Fits
# Open the PDF device
pdf("spline_fits_with_rgb_and_metrics.pdf", width = 12, height = 8)
# Set up the plotting area to show two plots per page
par(mfrow = c(1, 2))
# Loop over each species and plot the original data and spline fits side by side
for (species in names(spline_results)) {
result <- spline_results[[species]]
# Get the correct RGB color
color <- rgb_colors[grep(species, names(rgb_colors))]
# Plot original data and natural spline fit
plot(tanData_filtered$wl, tanData_filtered[[species]], main = paste("Natural Spline Fit for", species),
xlab = "Wavelength (nm)", ylab = "Reflectance", type = "o", col = color)
lines(tanData_filtered$wl, predict(result$natural_spline_fit), col = "red", lwd = 2)
mtext(paste("Adj R2:", round(result$natural_adj_r2, 3), "AIC:", round(result$natural_aic, 2)),
side = 3, line = 0.5, cex = 0.7, col = "blue")
# Plot original data and cubic spline fit
plot(tanData_filtered$wl, tanData_filtered[[species]], main = paste("Cubic Spline Fit for", species),
xlab = "Wavelength (nm)", ylab = "Reflectance", type = "o", col = color)
lines(tanData_filtered$wl, predict(result$cubic_spline_fit), col = "green", lwd = 2)
mtext(paste("Adj R2:", round(result$cubic_adj_r2, 3), "AIC:", round(result$cubic_aic, 2)),
side = 3, line = 0.5, cex = 0.7, col = "blue")
}
# Close the PDF device
dev.off()## quartz_off_screen
## 2Step 6 Print a Summary Table of Spline Fits
# Create the spline summary table with coefficients and goodness-of-fit metrics
spline_table <- data.frame(
Species = names(spline_results),
Natural_Spline_Coefficients = sapply(spline_results, function(x) paste(x$natural_coef, collapse = ", ")),
Cubic_Spline_Coefficients = sapply(spline_results, function(x) paste(x$cubic_coef, collapse = ", ")),
Natural_Adj_R2 = sapply(spline_results, function(x) round(x$natural_adj_r2, 3)),
Cubic_Adj_R2 = sapply(spline_results, function(x) round(x$cubic_adj_r2, 3)),
Natural_AIC = sapply(spline_results, function(x) round(x$natural_aic, 2)),
Cubic_AIC = sapply(spline_results, function(x) round(x$cubic_aic, 2))
)
# Remove row names
rownames(spline_table) <- NULL
# Display the table in the document without row names
kable(spline_table, caption = "Summary of Spline Fits for Each Species (Including Goodness-of-Fit)") %>%
kable_styling(bootstrap_options = "striped", full_width = F)| Species | Natural_Spline_Coefficients | Cubic_Spline_Coefficients | Natural_Adj_R2 | Cubic_Adj_R2 | Natural_AIC | Cubic_AIC |
|---|---|---|---|---|---|---|
| TanCyo | 5.49, 1.85, 6.47, 0.09, 2.03, 0.44 | 9, -0.07, 0, 0 | 0.920 | 0.717 | 434.02 | 940.85 |
| TanLab | 8.82, -2.1, 8.11, 1.57, -3.02, 5.53 | 124.37, -0.81, 0, 0 | 0.971 | 0.769 | 621.19 | 1447.47 |
| TanRue | 3.55, 0.23, 29.51, 5.17, 10.86, 2.65 | 177.07, -1.28, 0, 0 | 0.970 | 0.770 | 1229.01 | 2045.23 |
| TanChi | -4.73, -26.22, 79.79, 28.04, 39.68, -0.95 | 902.87, -6.16, 0.01, 0 | 0.867 | 0.659 | 2779.55 | 3156.06 |
| TanVel | 4.67, 3.79, -0.58, 0.53, -2.32, 2.04 | -25.97, 0.17, 0, 0 | 0.854 | 0.317 | 419.45 | 1034.89 |
| TanCal | 5.14, 20.49, 66.62, 86.72, 56.18, 37.57 | 812.16, -5.72, 0.01, 0 | 0.995 | 0.984 | 1739.27 | 2159.82 |
| TanMex | 2.83, 1.49, 1.72, 1.33, 3.24, 1.8 | -13.24, 0.1, 0, 0 | 0.851 | 0.848 | -34.11 | -29.58 |
| TanIno | 5.31, 4.18, 4.8, 5.56, 8.88, 6.11 | -17.72, 0.13, 0, 0 | 0.983 | 0.983 | 82.00 | 89.03 |
| TanSel | 2.54, -4.87, 117.83, 13.75, 31.09, -7.3 | 783.68, -5.71, 0.01, 0 | 0.925 | 0.715 | 2730.48 | 3265.24 |
| TanFas | 1.31, 2.88, 112.54, 10.91, 28.07, 0.89 | 634.65, -4.76, 0.01, 0 | 0.956 | 0.743 | 2471.70 | 3176.45 |
| TanCyc | 37.56, 12.07, -18.39, -27.75, -32.97, -24.42 | -323.85, 2.31, 0, 0 | 0.969 | 0.911 | 1788.56 | 2206.97 |
| TanCyv | 0.77, -4.97, 21.78, 29.97, 37.82, 27.44 | 274.9, -1.78, 0, 0 | 0.985 | 0.943 | 1483.02 | 2024.56 |
| TanDes | -1.07, -19.37, 66.97, 2.65, 28.12, -0.49 | 509.47, -3.53, 0.01, 0 | 0.902 | 0.550 | 2331.92 | 2942.65 |
| TanChr | 1.29, 3.58, 3.84, 8.51, 6.95, 3.27 | 34.67, -0.25, 0, 0 | 0.978 | 0.931 | 260.33 | 713.45 |
| TanXac | 3.57, -0.7, 2.63, 59.27, 56.46, 50.7 | 301.81, -1.81, 0, 0 | 0.995 | 0.939 | 1482.10 | 2490.76 |
| TanSch | -0.57, -9.88, 39.78, 44.33, 57.75, 42.84 | 471.6, -3.1, 0.01, 0 | 0.980 | 0.935 | 1968.85 | 2440.76 |
| TanJoh | 1.93, -14.12, 49.87, 12.8, 25.72, 3.14 | 474.41, -3.23, 0.01, 0 | 0.920 | 0.678 | 2112.46 | 2668.77 |
| TanArt | 0.23, -3.98, 13.2, 36.54, 41.89, 35.77 | 250.37, -1.57, 0, 0 | 0.993 | 0.960 | 1345.99 | 2017.60 |
| TanIct | -0.48, -9.35, 34, 33.51, 43.72, 27.5 | 410.4, -2.71, 0.01, 0 | 0.971 | 0.904 | 1889.41 | 2366.07 |
| TanFlo | -1.53, -24.02, 70.46, 59.26, 72.43, 41.1 | 927.96, -6.16, 0.01, 0 | 0.960 | 0.871 | 2512.40 | 2981.26 |
| TanPar | 3.34, -7.99, 33.57, 53.71, 53.47, 42.64 | 559.29, -3.66, 0.01, 0 | 0.993 | 0.959 | 1595.16 | 2304.96 |
| TanLav | 2.05, 0.16, -0.32, 12.54, 24.61, 32.5 | -58.09, 0.49, 0, 0 | 0.999 | 0.996 | 50.66 | 697.53 |
| TanGyr | 4.04, 0.37, -0.54, 10.6, 19.7, 25.76 | -40.58, 0.37, 0, 0 | 0.997 | 0.992 | 416.94 | 824.40 |
| TanVas | 4.47, 27.81, 20.65, 17.18, 10.92, 22.3 | -98.7, 0.4, 0, 0 | 0.979 | 0.817 | 1375.32 | 2239.01 |
| TanNiv | 6.08, 9.48, 108.33, 41.1, 29.4, 15.17 | 875.62, -6.37, 0.01, 0 | 0.989 | 0.883 | 2025.91 | 2956.90 |
| TanFuc | 4.26, 3.87, 44.72, 14.14, 15.7, 12.45 | 297.85, -2.18, 0.01, 0 | 0.995 | 0.861 | 939.51 | 2282.45 |
| TanDow | 9.29, -4.15, 5.65, 8.11, 4.92, 14.15 | 149.15, -0.92, 0, 0 | 0.987 | 0.972 | 731.23 | 1026.97 |
| TanRuu | 5.09, 1.49, 3.27, 5.51, 6.82, 8.33 | 17.61, -0.09, 0, 0 | 0.948 | 0.947 | 739.21 | 743.41 |
| TanPun | 3.82, 7.76, 35.69, 28.25, 22.32, 5.18 | 321.87, -2.3, 0.01, 0 | 0.955 | 0.919 | 1806.05 | 2040.05 |
| TanXag | 5.57, -29.14, 57.48, 4.05, 24.95, -8.26 | 617.88, -4.1, 0.01, 0 | 0.941 | 0.518 | 2071.70 | 2908.52 |
| TanGut | 2.26, -8.05, 35.49, 12.61, 23.07, 3.8 | 319.67, -2.17, 0, 0 | 0.956 | 0.737 | 1607.22 | 2318.38 |
| TanVar | 1.57, -13.14, 39.81, 4.7, 21.3, -2.47 | 331.05, -2.24, 0, 0 | 0.940 | 0.561 | 1717.09 | 2512.59 |
| TanPal | 15.54, 3.54, 5.36, 0, 14.83, -2.49 | -102.78, 0.78, 0, 0 | 0.863 | 0.605 | 891.24 | 1315.52 |
| TanVit | 2.25, 2.13, 5.97, 38.59, 39.52, 38.61 | 168.59, -1.03, 0, 0 | 0.999 | 0.969 | 661.58 | 1912.92 |
| TanCuc | 5.49, 0.1, 0.77, 2.67, 5.01, 6.22 | -0.38, 0.05, 0, 0 | 0.990 | 0.990 | -245.67 | -237.98 |
| TanCay | 3.33, 0.12, 8.37, 27.98, 28.83, 29.58 | 165.23, -1.03, 0, 0 | 0.999 | 0.984 | 107.89 | 1436.42 |
| TanPre | 9.87, -6.53, 1.96, 43.72, 25.99, 27.87 | 430.45, -2.69, 0.01, 0 | 0.995 | 0.938 | 1169.47 | 2211.25 |
| TanNic | 26.26, 42.67, 23.41, 0.32, 46.88, -15.12 | -680.49, 4.43, -0.01, 0 | 0.994 | 0.989 | 1051.49 | 1327.63 |
| TanCya | 10.54, 0.71, 81.85, -10.46, 53.36, -36.27 | 47.12, -0.51, 0, 0 | 0.977 | 0.661 | 1882.47 | 2963.59 |
| TanLar | 9.69, 14.72, 14.22, 74.71, 67.17, 80.11 | 288.47, -1.82, 0, 0 | 0.989 | 0.975 | 2016.22 | 2358.96 |
| TanCyp | 5, 0.16, 0.48, 1.2, 1.1, 2.38 | 10.89, -0.04, 0, 0 | 0.904 | 0.893 | -49.47 | -6.38 |
| TanVir | 6.5, -0.97, -0.72, 0.13, -1.12, 1.66 | 22.06, -0.1, 0, 0 | 0.868 | 0.824 | -105.50 | 7.12 |
| TanHei | 3.85, 0.96, 2.17, 1.88, 3.42, 2.48 | 0.76, 0.01, 0, 0 | 0.767 | 0.762 | 500.96 | 507.78 |
| TanArg | 2.59, 0.24, 0.56, 0.82, 0.87, 0.74 | 6.28, -0.03, 0, 0 | 0.616 | 0.616 | -28.33 | -30.27 |
| ThrCyp | 10.82, -3.61, 17.18, -0.92, 9.05, -2.02 | 92.98, -0.58, 0, 0 | 0.955 | 0.485 | 812.06 | 1786.58 |
| ThrAbb | 11.25, 7.92, 4.52, -11.19, 46.2, -25.53 | -503.23, 3.51, -0.01, 0 | 0.980 | 0.780 | 1535.84 | 2493.80 |
| ThrPal | 12.04, -16.38, 18.35, -8.19, 43.63, -15.23 | -69.59, 0.72, 0, 0 | 0.950 | 0.238 | 1581.24 | 2670.14 |
| ThrOrn | 26.98, 1.08, -10.46, -15.61, 4.52, -21.38 | -258.72, 1.97, 0, 0 | 0.982 | 0.930 | 1290.36 | 1821.61 |
| ThrSay | 7.51, 5.2, 18.59, 11.41, 24.03, 1.39 | 28.85, -0.2, 0, 0 | 0.976 | 0.783 | 732.42 | 1614.16 |
| ThrEpi | 17.47, 7.34, 36.49, -1.53, 31.31, -2.88 | -83.1, 0.53, 0, 0 | 0.969 | 0.664 | 1193.19 | 2149.73 |
| TanRuf | 2.9, 2.46, 8.23, 15.08, 17.93, 19.02 | 63.51, -0.41, 0, 0 | 0.999 | 0.997 | -100.03 | 243.82 |
# Save the table as a CSV
write.csv(spline_table, "spline_summary_table.csv", row.names = FALSE)
# Fixing the PDF export by increasing the margins and adjusting the layout
pdf("spline_summary_table.pdf", width = 11, height = 8.5)
gridExtra::grid.table(spline_table, rows = NULL)
dev.off()## quartz_off_screen
## 2Step 7 Save data objects for beast run
# Step 7: Save the Tanager Tree and Data as an R Data Object
# Save the tree, filtered spectral data, and spline results into an RData file
save(tanagerTree, tanData_filtered, spline_results, file = "tanager_analysis_results.RData")
# Inform the user that the data has been saved
cat("The Tanager tree, data, and spline results have been saved to 'tanager_analysis_results.RData'.\n")## The Tanager tree, data, and spline results have been saved to 'tanager_analysis_results.RData'.Step 8 Getting things ready for BEAST
# Export spline coefficients as a CSV file for BEAST
spline_traits <- data.frame(Species = names(spline_results))
# Add the coefficients as traits
for (i in 1:length(spline_results)) {
spline_traits[i, "Natural_Spline_1"] = spline_results[[i]]$natural_coef[2]
spline_traits[i, "Natural_Spline_2"] = spline_results[[i]]$natural_coef[3]
spline_traits[i, "Natural_Spline_3"] = spline_results[[i]]$natural_coef[4]
spline_traits[i, "Natural_Spline_4"] = spline_results[[i]]$natural_coef[5]
spline_traits[i, "Cubic_Spline_1"] = spline_results[[i]]$cubic_coef[2]
spline_traits[i, "Cubic_Spline_2"] = spline_results[[i]]$cubic_coef[3]
spline_traits[i, "Cubic_Spline_3"] = spline_results[[i]]$cubic_coef[4]
}
# Save the spline traits to a CSV file
write.csv(spline_traits, "spline_coefficients_for_BEAST.csv", row.names = FALSE)Step 9 BEAST xml
not sure yet how to do this but here is the general appropach….
Load Tree and Data: Load your phylogenetic tree (e.g., tanagerTree) and the CSV file containing the traits. Define the Traits: In the XML file, define each spline coefficient as a separate trait. If you want to estimate a single evolutionary rate for all coefficients, you can define them as independent traits. Model Specification: For each trait, specify a Brownian motion model (bm). If you want to treat the coefficients independently but estimate the same evolutionary rate for all, use a shared rate parameter across the traits. Rate Estimation: Define a hyperparameter for the rate of evolution and link it to all the traits. Alternatively, run separate BEAST analyses for each trait if you want to estimate the rates independently.
<trait id="natural_spline_1" spec="RealParameter" dimension="N">
<statefile name="state" startstate="true"/>
</trait>
<distribution id="prior" spec="util.CompoundDistribution">
<distribution spec="beast.math.distributions.MRCAPrior" id="TMRCA">
<tree topologies="$tanagerTree" spec="TreeLikelihood">
<distribution spec="beast.evolution.likelihood.MG94Likelihood" id="treeModel">
<data id="natural_spline_1.data" spec="beast.evolution.alignment.Alignment" dataType="continuous">
<sequence taxon="species1" totalcount="4" value="coeff1, coeff2, coeff3, coeff4"/>
<sequence taxon="species2" totalcount="4" value="coeff1, coeff2, coeff3, coeff4"/>
</data>
</distribution>
</tree>
</distribution>
</distribution>
<run id="mcmc" spec="MCMC">
<state>
<parameter id="rate" name="stateNode" dimension="1">0.1</parameter>
</state>
<posterior>
<distribution id="posterior" spec="util.CompoundDistribution">
<prior id="prior" spec="beast.math.distributions.MRCAPrior" tree="@treeModel">
<normalPrior mean="0" stdev="1" name="distr"/>
</prior>
<likelihood id="likelihood" spec="util.CompoundDistribution">
<distribution id="treeLikelihood" spec="TreeLikelihood">
<data id="alignment" spec="Alignment" dataType="continuous"/>
<siteModel id="SiteModel" spec="SiteModel" gammaCategoryCount="4">
<substitutionModel id="F81" spec="F81"/>
</siteModel>
</distribution>
</likelihood>
</distribution>
</posterior>
</run>You will need to replicate this setup for each coefficient, ensuring that they either share or have independent rate parameters depending on your choice.
- Run BEAST
Load the XML into BEAST: Once the XML is set up, load it into BEAST. Run the Analysis: Start the analysis to estimate the evolutionary rates. 4. Post-Analysis
After the BEAST run completes, you can analyze the results using Tracer to check the estimated rates and other parameters.